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Abstract 

In this paper, we model the evolution and self-assembly of randomly oriented 
carbon nanotubes (CNTs), grown on a metallic substrate in the form of a thin film 
for field emission under diode configuration. Despite high output, the current in 
such a thin film device often decays drastically. The present paper is focused on 
understanding this problem. A systematic, multiphysics based modelling approach 
is proposed. First, a nucleation coupled model for degradation of the CNT thin film 
is derived, where the CNTs are assumed to decay by fragmentation and formation of 
clusters. The random orientation of the CNTs and the electromechanical interaction 
are then modeled to explain the self-assembly. The degraded state of the CNTs and 
the electromechanical force are employed to update the orientation of the CNTs. 
Field emission current at the device scale is finally obtained by using the Fowler- 
Nordheim equation and integration over the computational cell surfaces on the anode 
side. The simulated results are in close agreement with the experimental results. 
Based on the developed model, numerical simulations aimed at understanding the 
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effects of various geometric parameters and their statistical features on the device 
current history are reported. 

Keywords: Field emission, carbon nanotube, degradation, electrodynamics, self- 
assembly. 



1 Introduction 



The conventional mechanism used for electron emission is thermionic in nature where 
electrons are emitted from hot cathodes (usually heated filaments). The advantage 
of these hot cathodes is that they work even in environments that contain a large 
number of gaseous molecules. However, thermionic cathodes in general have slow 
response time and they consume high power. These cathodes have limited lifetime 
due to mechanical wear. In addition, the thermionic electrons have random spatial 
distribution. As a result, fine focusing of electron beam is very difficult. This 
adversely affects the performance of the devices such as X-ray tubes. An alternative 
mechanism to extract electrons is field emission, in which electrons near the Fermi 
level tunnel through the energy barrier and escape to the vacuum under the influence 
of a sufficiently high external electric field. The field emission cathodes have faster 
response time, consume less power and have longer life compared to thermionic 
cathodes. However, field emission cathodes require ultra-high vacuum as they are 
highly reactive to gaseous molecules during the field emission. 

The key to the high performance of a field emission device is the behavior of 
its cathode. In the past, the performance of cathode materials such as spindt-type 
emitters and nanostructured diamonds for field emission was studied by Spindt et 
al. 1 , Gotoh et al. 2 , and Zhu 3 . However, the spindt type emitters suffer from 
high manufacturing cost and limited lifetime. Their failure is often caused by ion 
bombardment from the residual gas species that blunt the emitter cones 2 . On 
the other hand, nanostructured diamonds are unstable at high current densities 3 . 
Carbon nanotube (CNT), which is an allotrope of carbon, has potential to be used 
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as cathode material in field emission devices. Since their discovery by Iijima in 1991 4 
, extensive research on CNTs has been conducted. Field emission from CNTs was 
first reported in 1995 by Rinzler et al. 5 , de Heer et al. 6 , and Chernozatonskii et 
al. 7 . Field emission from CNTs has been studied extensively since then. Currently, 
with significant improvement in processing technique, CNTs are among the best field 
emitters. Their applications in field emission devices, such as field emission displays, 
gas discharge tubes, nanolithography systems, electron microscopes, lamps, and X- 
ray tube sources have been successfully demonstrated 8-9 . The need for highly 
controlled application of CNTs in X-ray devices is one of the main reasons for the 
present study. The remarkable field emission properties of CNTs are attributed 
to their geometry, high thermal conductivity, and chemical stability. Studies have 
reported that CNT sources have a high reduced brightness and their energy spread 
values are comparable to conventional field emitters and thermionic emitters 10 . 

The physics of field emission from metallic surfaces is well understood. The 
current density (J) due to field emission from a metallic surface is usually obtained 
by using the Fowler-Nordheim (FN) equation 11 



where E is the electric field, $ is the work function of the cathode material, and B 
and C are constants. The device under consideration in this paper is a X-ray source 
where a thin film of CNTs acts as the electron emitting surface (cathode). Under the 
influence of sufficiently high voltage at ultra high vacuum, the electrons are extracted 
from the CNTs and hit the heavy metal target (anode) to produce X-rays. However, 
in the case of a CNT thin film acting as cathode, the surface of the cathode is not 
smooth (like the metal emitters). In this case, the cathode consists of hollow tubes 
grown on a substrate. Also, some amount of carbon clusters may be present within 
the CNT-based film. An added complexity is that there is realignment of individual 
CNTs due to electrodynamic interaction between the neighbouring CNTs during 
field emission. At present, there is no adequate mathematical models to address 
these issues. Therefore, the development of an appropriate mathematical modeling 
approach is necessary to understand the behavior of CNT thin film field emitters. 




(1) 
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1.1 Role of various physical processes in the degradation of 
CNT field emitter 

Several studies have reported experimental observations in favour of considerable 
degradation and failure of CNT cathodes. These studies can be divided into two 
categories: (i) studies related to degradation of single nanotube emitters 12-17 and 
(ii) studies related to degradation of CNT thin films 18-23 . Dean et al. 13 found 
gradual decrease of field emission of single walled carbon nanotubes (SWNTs) due 
to "evaporation" when large field emitted current (300nA to 1\xA) was extracted. It 
was observed by Lim et al. 23 that CNTs are susceptible to damage by exposure to 
gases such as oxygen and nitrogen during field emission. Wei et al. 14 observed that 
after field emission over 30 minutes at field emission current between 50 and 120nA, 
the length of CNTs reduced by 10%. Wang et al. 15 observed two types of structural 
damage as the voltage was increased: a piece-by-piece and segment-by-segment 
splitting of the nanotubes, and a layer-by-layer stripping process. Occasional spikes 
in the current- voltage curves were observed by Chung et al. 16 when the voltage 
was increased. Avouris et al. 17 found that the CNTs break down when subjected 
to high bias over a long period of time. Usually, the breakdown process involves 
stepwise increases in the resistance. In the experiments performed by the present 
authors, peeling of the film from the substrate was observed at high bias. Some of 
the physics pertinent to these effects is known but the overall phenomenon governing 
such a complex system is difficult to explain and quantify and it requires further 
investigation. 

There are several causes of CNT failures: 

(i) In case of multi- walled carbon nanotubes (MWNTs), the CNTs undergo layer- 
by-layer stripping during field emission 15 . The complete removal of the shells 
are most likely the reason for the variation in the current voltage curves 16 ; 

(ii) At high emitted currents, CNTs are resistively heated. Thermal effect can 
sublime a CNT causing cathode-initiated vacuum breakdown 24 . Also, in case 
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of thin films grown using chemical vapor deposition (CVD), fewer catalytic 
metals such as nickel, cobalt, and iron are observed as impurities in CNT thin 
films. These metal particles melt and evaporate by high emission currents, 
and abruptly surge the emission current. This results in vacuum breakdown 
followed by the failure of the CNT film 23 ; 

(iii) Gas exposure induces chemisorption and physisorption of gas molecules on the 
surface of CNTs. In the low-voltage regime, the gas adsorbates remain on the 
surface of the emitters. On the other hand, in the high-voltage regime, large 
emission currents resistively anneal the tips, and the strong electric field on the 
locally heated tips promotes the desorption of gas adsorbates from the tip sur- 
face. Adsorption of materials with high electronegativity hinders the electron 
emission by intensifying the local potential barriers. Surface morphology can 
be changed by an erosion of the cap of the CNT as the gases desorb reactively 
from the surface of the CNTs 25 ; 

(iv) CVD-grown CNTs tend to show more defects in the wall as their radius in- 
creases. Possibly, there are rearrangements of atomic structures (for example, 
vacancy migration) resulting in the reduction of length of CNTs 14 . In addition, 
the presence of defects may act as a centre for nucleation for voltage-induced 
oxidation, resulting in electrical breakdown 16 ; 

(v) As the CNTs grow perpendicular to the substrate, the contact area of CNTs 
with the substrate is very small. This is a weak point in CNT films grown on 
planar substrates, and CNTs may fail due to tension under the applied fields 20 
. Small nanotube diameters and lengths are an advantage from the stability 
point of view. 

Although the degradation and failure of single nanotube emitters can be either 
abrupt or gradual, the degradation and failure of a thin film emitter with CNT clus- 
ter is mostly gradual. The gradual degradation occurs either during initial current- 
voltage measurement 21 (at a fast time scale) or during measurements at constant 
applied voltage over a long period of time 22 (at a slow time scale). Nevertheless, it 



5 



can be concluded that the gradual degradation of thin films occurs due to the failure 
of individual emitters. 

Till date, several studies have reported experimental observations on CNT thin 
films 26 . However, from mathematical, computational and design view points, the 
models and characterization methods are available only for vertically aligned CNTs 
grown on the patterned surface 27-28 . In a CNT film, the array of CNTs may ideally 
be aligned vertically However, in this case it is desired that the individual CNTs 
be evenly separated in such a way that their spacing is greater than their height to 
minimize the screening effect 29 . If the screening effect is minimized following the 
above argument, then the emission properties as well as the lifetime of the cathodes 
are adversely affected due to the significant reduction in density of CNTs. For the 
cathodes with randomly oriented CNTs, the field emission current is produced by 
two types of sources: (i) small fraction of CNTs that point toward the anode and 
(ii) oriented and curved CNTs subjected to electromechanical forces causing reori- 
entation. As often inferred (see e.g., ref. 29 ), the advantage of the cathodes with 
randomly oriented CNTs is that always a large number of CNTs take part in the field 
emission, which is unlikely in the case of cathodes with uniformly aligned CNTs. 
Such a thin film of randomly oriented CNTs will be considered in the present study. 
From the modeling point of view, its analysis becomes much more challenging. Al- 
though some preliminary works have been reported (see e.g., refs. 30-31 ) ; neither a 
detailed model nor a subsequent characterization method are available that would 
allow to describe the array of CNTs that may undergo complex dynamics during 
the process of charge transport. In the detailed model, the effects of degradation 
and fragmentation of CNTs during field emission need to be considered. However, 
in the majority of analytical and design studies, the usual practice is to employ 
the classical Fowler-Nordheim equation 11 to determine the field emission from the 
metallic surface, with correction factors to deal with the CNT tip geometry. Ideally, 
one has to tune such an empirical approach to specific materials and methods used 
(e.g. CNT geometry, method of preparation, CNT density, diode configuration, 
range of applied voltage, etc.). Also, in order to account for the oriented CNTs and 
interaction between themselves, it is necessary to consider the space charge and the 
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electromechanical forces. By taking into account the evolution of the CNTs, a mod- 
eling approach is developed in this paper. In order to determine phenomenologically 
the concentration of carbon clusters due to degradation of CNTs, we introduce a 
homogeneous nucleation rate. This rate is coupled to a moment model for the evo- 
lution. The moment model is incorporated in a spatially discrete sense, that is by 
introducing volume elements or cells to physically represent the CNT thin film. Elec- 
tromechanical forces acting on the CNTs are estimated in time-incremental manner. 
The oriented state of CNTs are updated using a mechanics based model. Finally, 
the current density is calculated by using the details regarding the CNT orientation 
angle and the effective electric field in the Fowler-Nordheim equation. 

The remainder of this paper is organized as follows: in Sec. 2, a model is pro- 
posed, which combines the nucleation coupled model for CNT degradation with the 
electromechanical forcing model. Section 3 illustrates the computational scheme. 
Numerical simulations and the comparison of the simulated current- volt age charac- 
teristics with experimental results are presented in Sec. 4. 

2 Model formulation 

The CNT thin film is idealized in our mathematical model by using the following 
simplifications. 

(i) CNTs are grown on a substrate to form a thin film. They are treated as 
aggregate while deriving the nucleation coupled model for degradation phe- 
nomenologically; 

(ii) The film is discretized into a number of representative volume element (cell), 
in which a number of CNTs can be in oriented forms along with an estimated 
amount of carbon clusters. This is schematically shown in Fig. [T] The car- 
bon clusters are assumed to be in the form of carbon chains and networks 
(monomers and polymers); 
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(iii) Each of the CNTs with hexagonal arrangement of carbon atoms (shown in 
Fig. [21(a)) are treated as effectively one-dimensional (ID) elastic members and 
discretized by nodes and segments along its axis as shown in Fig. [2^b). Defor- 
mation of this ID representation in the slow time scale defines the orientations 
of the segments within the cell. A deformation in the fast time scale (due to 
electron flow) defines the fluctuation of the sheet of carbon atoms in the CNTs 
and hence the resulting state of atomic arrangements. The latter aspect is ex- 
cluded from the present modeling and numerical simulations, however they 
will be discussed within a quantum-hydrodynamic framework in a forthcom- 
ing article. 

2.1 Nucleation coupled model for degradation of CNTs 

Let Nt be the total number of carbon atoms (in CNTs and in cluster form) in a 
cell (see Fig. [I]). The volume of a cell is given by V ce \\ = A Ad, where A A is the cell 
surface interfacing the anode and d is distance between the inner surfaces of cathode 
substrate and the anode. Let N be the number of CNTs in the cell, and ./Vcnt be 
the total number of carbon atoms present in the CNTs. We assume that during 
field emission some CNTs are decomposed and form clusters. Such degradation and 
fragmentation of CNTs can be treated as the reverse process of CVD or a similar 
growth process used for producing the CNTs on a substrate. Hence, 



where -/Venter is the total number of carbon atoms in the clusters in a cell at time t 
and is given by 



where n\ is the concentration of carbon cluster in the cell. By combining Eqs. rt2 



N T = NN CNT + N clustcI , 



(2) 




(3) 



and one has 
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The number of carbon atom in a CNT is proportional to its length. Let the length 
of a CNT be a function of time, denoted as L(t). Therefore, one can write 

Ncnt = N ving L(t) , (5) 

where N ring is the number of carbon atoms per unit length of a CNT and can be 
determined from the geometry of the hexagonal arrangement of carbon atoms in the 
CNT. By combining Eqs. Q and ([5]), one can write 

1 



N 



N ving L(t) 



t 

N T -V ccll I dn x {t) . (6) 
o 



In order to determine n\(t) phenomenologically, we need to know the nature of 
evolution of the aggregate in the cell. From the physical point of view, one may 
expect the rate of formation of the carbon clusters from CNTs to be a function 
of thermodynamic quantities, such as temperature (T), the relative distances (ry) 
between the carbon atoms in the CNTs, the relative distances between the clusters 
and a set of parameters (p*) describing the critical cluster geometry. The relative 
distance between carbon atoms in CNTs is a function of the electromechanical 



forces. Modeling of this effect is discussed in Sec. |2.2[ On the other hand, the relative 
distances between the clusters influence in homogenizing the thermodynamic energy, 
that is, the decreasing distances between the clusters (hence increasing densities of 
clusters) slow down the rate of degradation and fragmentation of CNTs and lead to 
a saturation in the concentration of clusters in a cell. Thus, one can write 

^=/(T, W *). (7) 

To proceed further, we introduce a nucleation coupled model 32-33 , which was origi- 
nally proposed to simulate aerosol formation. Here we modify this model according 
to the present problem which is opposite to the process of growth of CNTs from 
the gaseous phase. With this model the relative distance function is replaced by a 
collision frequency function (fly) describing the frequency of collision between the 
z-mers and j-mers, with 



and the set of parameters describing the critical cluster geometry by 



p* = { Vj Sj g* d* p } , (9) 

where Vj is the j-mer volume, Sj is the surface area of j-mer, g* is the normalized 
critical cluster size, d* is the critical cluster diameter, k is the Boltzmann constant, T 
is the temperature and p p is the particle mass density. The detailed form of Eq. ([7| 
is given by four nonlinear ordinary differential equations: 

dN kin 



dt 



4in , (10) 



dS JkmSg* 
dM x 



dt 



Jkin< + (S - l^iiVkin , (12) 



dA n _ J kin Sg* 2 / 3 Sl 2-kB 1 S(S-1)M 1 

dt Til 7l\ 

where jVyi n is the kinetic normalization constant, Jkm is the kinetic nucleation rate, 
S is the saturation ratio, A n is the total surface area of the carbon cluster and Mi 
is the moment of cluster size distribution. The quantities involved are expressed as 

3= — , Mi = I (n(d p ,t)d p )d(d p ) , (14) 

Tic I V / 



dp 



Man = -i exp(6) , J kin = ^-f J — exp - — — , (15) 



S ^ ' ' Km 125 V 2tt " V 27(lnS) 2 



/2 \ 3 I!k 4av! kT 

where n s is the equilibrium saturation concentration of carbon cluster, <i™ ax is the 
maximum diameter of the clusters, n(d p , t) is the cluster size distribution function, 
dp is the cluster diameter, rrij is the mass of j-mer, is the dimensionless surface 
tension given by 

a is the surface tension. In this paper, we have considered 2 = 1 and j = 1 for numer- 



ical simulations, that is, only monomer type clusters are considered. In Eqs. (10)- 



(13), the variables are rti(i), S(t), M\(t) and A n (t), and all other quantities are 
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assumed constant over time. In the expression for moment Mi(t) in Eq. (14), the 
cluster size distribution in the cell is assumed to be Gaussian, however, random 



distribution can be incorporated. We solve Eqs. ( 10 )-( 13 ) using a finite difference 
scheme as discussed in Sec. [3} Finally, the number of CNTs in the cell at a given time 
is obtained with the help of Eq. (|6]), where the reduced length Lit) is determined 
using geometric properties of the individual CNTs as formulated next. 



2.2 Effect of CNT geometry and orientation 



It has been discussed in Sec. |1.1| that the geometry and orientation of the tip of the 
CNTs are important factors in the overall field emission performance of the film and 
must be considered in the model. 

As an initial condition, let L(0) = h at t = 0, and let ho be the average height of 
the CNT region as shown in Fig. [TJ This average height ho is approximately equal 
to the height of the CNTs that are aligned vertically. If Ah is the decrease in the 
length of a CNT (aligned vertically or oriented as a segment) over a time interval 
At due to degradation and fragmentation, and if dt is the diameter of the CNT, 
then the surface area of the CNT decreased is nd t Ah. By using the geometry of the 
CNT, the decreased surface area can be expressed as 

r i 1/2 

ird t Ah = V ce nni(t) s(s - ai)(s - a 2 )(s - a 3 ) , (18) 



where V ce n is the volume of the cell as introduced in Sec. |2.1[ a±, a 2 , a 3 are the lattice 
constants, and s — \{a>i + a 2 + a 3 ) (see Fig. |2j(a)). The chiral vector for the CNT is 
expressed as 

C h = nai + ma 2 , (19) 

where n and m are integers (n > \m\ > 0) and the pair (n,m) defines the chirality 
of the CNT. The following properties hold: d\.di = af, a 2 .a 2 = a|, and 2d\.a 2 = 
a\ + a\ — a\. With the help of these properties the circumference and the diameter 
of the CNT can be expressed as, respectively 34 , 



C h\ = \l n 2 a\ + m 2 a\ + nm{a\ + a| — a|) , d t = - — — , (20) 

7T 
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Let us now introduce the rate of degradation of the CNT or simply the burning rate 
as fbum = hn 

At- 

limit, one has 



as fbum = ^ m Ah /At. By dividing both side of Eq. (18) by At and by applying 



ndtVbum = V c 



ceir 



dni{t) 
dt 



1/2 

s(s - a 1 )(s - a 2 )(s - a 3 ) , (21) 



By combining Eqs. (20) and ( |21[ ), the burning rate is finally obtained as 

dn\{t) r s(s — ai)(s — a 2 )(s — a 3 ) i 1 / 2 



^ burn 



cell" 



n' z af + m?a\ + nm(a( + a| — a| 



(22) 



In Fig. [3] we show a schematic drawing of the CNTs almost vertically aligned, 
that is along the direction of the electric field E(x,y). This electric field E(x,y) is 
assumed to be due to the applied bias voltage. However, there will be an additional 
but small amount of electric field due to several localized phenomena (e.g., electron 
flow in curved CNTs, field emission from the CNT tip etc.). Effectively, we assume 
that the distribution of the field parallel to z-axis is of periodic nature (as shown 
in Fig. [3| when the CNT tips are vertically oriented. Only a cross-sectional view in 
the xz plane is shown in Fig. [3] because only an array of CNTs across x-direction 
will be considered in the model for simplicity. Thus, in this paper, we shall restrict 
our attention to a two-dimensional problem, and out-of-plane motion of the CNTs 
will not be incorporated in the model. 

To determine the effective electric field at the tip of a CNT oriented at an angle 
9 as shown in Fig. |3j we need to know the tip coordinate with respect to the cell 
coordinate system. If it is assumed that a CNT tip was almost vertically aligned at 
t — (as it is the desired configuration for the ideal field emission cathode), then 
its present height is L(t) = h — Vb um t an d the present distance between the tip 
and the anode is d g = d — L(t) = d — h + fburn£- We assume that the tip electric 
field has a ^-dependence of the form E L(t)/d g , where E = V/d and V is the 
applied bias voltage. Also, let (x, y) be the deflection of the tip with respect to its 
original location and the spacing between the two neighboring CNTs at the cathode 
substrate is 2R. Then the electric field at the deflected tip can be approximated as 
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where 9 C is a critical angle to be set during numerical calculations along with the 
condition: E z > = when 9(t) > 9 C . This is consistent with the fact that those CNTs 
which are low lying on the substrate do not contribute to the field emission. The 
electric field at the individual CNT tip derived here is defined in the local coordinate 
system (X', Z') as shown in Fig. [3] The components of the electric field in the cell 
coordinate system (X, Y, Z) is given by the following transformation: 



" E z ' 




E x 




. Ey . 











m z 




m z n z 






-hn z 


-lz 









E z i 













(24) 



where n z , l z , m z are the direction cosines. According to the cell coordinate system 
in Figs. [I] and [3j n z = cos0(t), l z = sm9(t), and m z = 0. Therefore, Eq. (24) can 
be rewritten as 



E z 




cos9(t) 


sm9(t) 




E z i 


E x 




y/1 ~ COS 2 6(t) 


- cos 0(t) 







Ey 







- cos 9(t) -1 








(25) 



By simplifying Eq. (25), we get 

E z = E z , cos 9 it) 



E x = E z , sin 9 (t) . 



(26) 



Note that the identical steps of this transformation also apply to a generally oriented 
(9 7^ 0) segment of CNT as idealized in Fig. |2^b). The electric field components E z 
and E x are later used for calculation of the electromechanical force acting on the 
CNTs. Since in this study we aim at estimating the current density at the anode due 



to the field emission from the CNT tips, we also use E z from Eq. (26) to compute 
the output current based on the Fowler- Nordheim equation (lip. 



2.3 Electromechanical forces 

For each CNT, the angle of orientation 9(f) is dependent on the electromechanical 
forces. Such dependence is geometrically nonlinear and it is not practical to solve 
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the problem exactly, especially in the present situation where a large number of 
CNTs are to be dealt with. However, it is possible to solve the problem in time- 
dependent manner with an incremental update scheme. In this section we derive 
the components of the electromechanical forces acting on a generally oriented CNT 
segment. The numerical solution scheme based on an incremental update scheme 
will be discussed in Sec. HI 

From the studies reported in published literature and based on the discussions 
made in Sec. it is reasonable to expect that the major contribution is due to (i) 
the Lorentz force under electron gas flow in CNTs (a hydrodynamic formalism), (ii) 
the electrostatic force (background charge in the cell), (iii) the van der Waals force 
against bending and shearing of MWNT and (iv) the ponderomotive force acting on 
the CNTs. 

2.3.1 Lorentz force 

It is known that the electrical conduction and related properties of CNTs depend on 
the mechanical deformation and the geometry of the CNT. In this paper we model 
the field emission behaviour of the CNT thin film by considering the time-dependent 
electromechanical effects, whereas the electronic properties and related effects are 
incorporated through the Fowler-Nordheim equation empirically. Electronic band- 
structure calculations are computationally prohibitive at this stage and at the same 
spatio-temporal scales considered for this study. However, a quantum-hydrodynamic 
formalism seems practical and such details will be dealt in a forthcoming article. 
Within the quantum-hydrodynamic formalism, one generally assumes the flow of 
electron gas along the cylindrical sheet of CNTs. The associated electron density 
distribution is related to the energy states along the length of the CNTs including the 
tip region. What is important for the present modeling is that the CNTs experience 
Lorentz force under the influence of the bias electric field as the electrons flow from 
the cathode substrate to the tip of a CNT. The Lorentz force is expressed as 

fi = e(n + n x )E « en E , (27) 
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where e is the electronic charge, n is the surface electron density corresponding to 
the Fermi level energy, hi is the electron density due to the deformation in the slow 
time scale, and phonon and electromagnetic wave coupling at the fast time scale, 
and E is the electric field. The surface electron density corresponding to the Fermi 
level energy is expressed as 35 

kT 

*> = as • ( 28 > 

where b is the interatomic distance and A is the overlap integral (~ 2eV for carbon). 
The quantity b can be related to the mechanical deformation of the ID segments (See 
Fig. [2]) and formulations reported by Xiao et al. 3e can be employed. For simplicity, 
the electron density fluctuation hi is neglected in this paper. Now, with the electric 
field components derived in Eq. [26j the components of the Lorentz force acting along 
z and x directions can now be written as, respectively, 

fi„ = %d t eh E z , f lx = nd t eh E x « . (29) 



2.3.2 Electrostatic force 

In order to calculate the electrostatic force, the interaction among two neighboring 
CNTs is considered. For such calculation, let us consider a segment dsi on a CNT 
(denoted 1) and another segment ds2 on its neighboring CNT (denoted 2). These 
are parts of the representative ID member idealized as shown in Fig. [2|b). The 
charges associated with these two segments can be expressed as 

qi = ehoirdt dsi , q 2 = en^ird^f 1 ds 2 , (30) 

where d[^ and d[ are diameters of two neighbouring CNTs (1) and (2). The 
electrostatic force on the segment dsi by the segment ds 2 is 

4vree rf 2 

where e is the effective permittivity of the aggregate of CNTs and carbon clusters, 
eo is the permittivity of free space, and r\% is the effective distance between the 
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centroids of dsi and ds 2 . The electrostatic force on the segment ds\ due to charge 
in the entire segment (s 2 ) of the neighboring CNT (see Fig. El) can be written as 



/ — (eh 7cd^ dsienQTxd!^\ ds 2 ■ 

J0 r l2 ^ ' 



fc = ~ A / " V ds * ■ ( 31 ) 



4vree j ri 2 

The electrostatic force per unit length on si due to s 2 is then 

1 n (7renofjM 2) 
4vree 7 rf 2 

The differential of the force d/ c acts along the line joining the centroids of the 
segments ds\ and ds 2 as shown in Fig. |4j Therefore, the components of the total 
electrostatic force per unit length of CNT (1) in X and Z directions can be written 
as, respectively, 

f 1 f s ' 2 (7ref?,o) 2 4 1) 4 2) 

fc x = / df c COS = / ~ COS ^ 

_ 1 ^ 7ren yd\ >d\> 

y 2 cos0 As 2 , (32) 



47ree ^ r 12 



/ ^ • a, 1 r 2 (^^o)^i 1} 4 2) . , , 

fc, = / dfc sin = / 5 sin ds 2 



4vree Jo rj 2 

= — g ^ sm0A, 2 , (33) 

where is the angle the force vector d/ c makes with the X-axis. For numerical 

computation of the above integrals, we compute the angle = 0(s^,s 2 ) and r 12 = 

r i2(<Si, 5 2 ) at each of the centroids of the segments between the nodes k + 1 and fc, 

where the length of the segments are assumed to be uniform and denoted as Asi 

for CNT (1) and As 2 for CNT (2). As shown in Fig. |1J the distance r% 2 between the 

centroids of the segments ds\ and ds 2 is obtained as 

r i 1/2 

r 12 = ^(di - l X2 + l Xl ) 2 + {1 Z1 - Z 22 ) 2 j , (34) 

where d\ is the spacing between the CNTs at the cathode substrate, l Xl and l X2 are 
the deflections along X-axis, and l Zl and l Z2 are the deflections along Z-axis. The 
angle of projection is expressed as 

-if '*1 ^2 



\di — l X2 + / Xl 
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The deflections 1 X1 , 1 Z1 , l X2 , and l Z2 are defined as, respectively, 



lx! — / dsi sin 9i = Asi sin 9{ 
Jo j 

/si 
ds\ cos 9i = Agj cos fl{ 

3 

/•S2 

i ffi2 = / <is 2 sin 9 2 = As2 sin 9 J 2 
Jo j 

ds 2 cos 9 2 = ^^As 2 cos^ . 



(36) 

(37) 
(38) 
(39) 



Note that the total electrostatic force on a particular CNT is to be obtained by 
summing up all the binary contributions within the cell, that is by summing up 



Eqs. (32) and (33) over the upper integer number of the quantity N — 1, where N 
is the number of CNTs in the cell as discussed in Sec. 12.11 



2.3.3 The van der Waals force 

Next, we consider the van der Waals effect. The van der Waals force plays important 
role not only in the interaction of the CNTs with the substrate, but also in the 
interaction between the walls of MWNTs and CNT bundles. Due to the overall 
effect of forces and flexibility of the CNTs (here assumed to be elastic ID members), 
the cylindrical symmetry of CNTs is destroyed, leading to their axial and radial 
deformations. The change in cylindrical symmetry may significantly affect the the 
properties of CNTs 37 " 38 . Here we estimate the van der Waals forces due to the 
interaction between two concentric walls of the MWCNTs. 

Let us assume that the lateral and the longitudinal displacements of a CNT be 
u x i and u z >, respectively. We use updated Lagrangian approach with local coordi- 
nate system for this description (similar to (X', Z') system shown in Fig. [3|, where 
the longitudinal axis coincides with Z' and the lateral axis coincides with X' . Such 
a description is consistent with the incremental procedure to update the CNT orien- 
tations in the cells as adopted in the computational scheme. Also, due to the large 
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length-to-diameter ratio (L(t)/d t ), let the kinematics of the CNTs, which are ideal- 
ized in this work as ID elastic members, be governed by that of an Euler-Bernoulli 
beam. Therefore, the kinematics can be written as 

r\ (m) 

u M = u H_ r (m)^_ , (40 ) 

where the superscript (m) indicates the mth wall of the MWNT with r^ m ' as its 
radius and u z 'q is the longitudinal displacement of the center of the cylindrical cross- 
section. Under tension, bending moment and lateral shear force, the elongation of 
one wall relative to its neighboring wall is 

Q„,( m + 1 ) P)oS m ^ £>/\ 

* (m) = (m+1) _ (m) = (m+D OU x> _ (m) OU x' ^ (Am+l) _ (m)\ °^' (at) 

dz' dz' ~ v ; ds ' v ; 

where we assume = u^ +1 ^ = A x > as the lateral displacement as some function 
of tensile force or compression buckling or pressure in the thin film device. The 
lateral shear stress (ris ) due to the van der Waals effect can now be written as 

A ( m ) 

rif =C VS -^, (42) 
where C vs is the van der Waals coefficient. Hence, the shear force per unit length 



can be obtained by integrating Eq. (42) over the individual wall circumferences and 
then by summing up for all the neighboring pair interactions, that is, 

f - = \], C ^ r ^ d ^ = \], C - { 2 J # 

=► Us = $>C4(r^) 2 - (r^) 2 ]^ 9 ^ • (43) 

m ' 

The components of van der Waals force in the cell coordinate system (X', Z') is then 
obtained as 

fvsz = f vs sm9(t) , f VSx = f vs cos 9(t) . (44) 



2.3.4 Ponderomotive force 

Ponderomotive force, which acts on free charges on the surface of CNTs, tends to 
straighten the bent CNTs under the influence of electric field in the Z-direction. 
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Furthermore, the ponderomotive forces induced by the applied electric field stretch 
every CNT 39 . We add this effect by assuming that the free charge at the tip region 
is subjected to Ponderomotive force, which is computed as 40 

f Pz = ^e E 2 AA cos 9(t) , f Px =0, (45) 

where AA is the surface area of the cell on the anode side, f Pz is the Z component 
of the Ponderomotive force and the X component f Px is assumed to be negligible. 



2.4 Modelling the reorientation of CNTs 

The net force components acting on the CNTs along Z and X directions can be 
expressed as, respectively, 

fz = J (fi z + fvJds + f Cz + f Pz , (46) 

U = J (fir ■ fvjds I I, I J),,.. (47) 
For numerical computation, at each time step the force components obtained using 



Eqs. (46) and (47) are employed to update the curved shape S'(x' + u x >,z' + u z /), 
where the displacements are approximated using simple beam mechanics solution: 

u z > « ^(Jt 1 - fi)(z' j+1 - , (48) 

« ^FJ-Ut 1 ~ fl) - , (49) 

where A is the effective cross-sectional area, A 2 is the area moment, E' is the 
modulus of elasticity for the CNT under consideration. The angle of orientation, 
9(t), of the corresponding segment of the CNT, that is between the node j + 1 and 
node j, is given by 

, w ^ (f) ^ tan -.( ^;+"i;:?-'i:t ). (50) 

V [z3 +1 + uV ) - [zi + ui) J 
ij I I ui, 



[V(9(t - AtY)] f , (51) 



ui " ui, 
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where T is the usual coordinate transformation matrix which maps the displace- 
ments (u X ',u z i) defined in the local (X',Z') coordinate system into the displace- 
ments (u x , u z ) defined in the cell coordinate system (X, Z). For this transformation, 
we employ the angle 9(t — At) obtained in the previous time step and for each node 
j = l,2,.... 



3 Computational scheme 

As already highlighted in the previous section, we model the CNTs as generally 
oriented ID elastic members. These ID members are represented by nodes and 
segments. With given initial distribution of the CNTs in the cell, we discretize 
the time into uniform steps tj + i — ti — At. The computational scheme involves 
three parts: (i) discretization of the nucleation coupled model for degradation of 



CNTs derived in Sec. 2.1, (ii) incremental update of the CNT geometry using the 
estimated electromechanical force and (iii) computation of the field emission current 
in the device. 



3.1 Discretization of the nucleation coupled model for degra- 
dation 



With the help of Eqs. (14)-(16) and by eliminating the kinetic nucleation rate N] 



kin; 



we first rewrite the simplified form of Eqs. (jlOp— (jl3p , which are given by, respectively, 

46 3 



S 



ni- 



dS 
~dt 



dn\ 

dt x dt 

2f3 lins QS 
' 81\/2^(ln sy 



dS pnn 2 s S 3 Q 



12 



cxp 



dMi Pnn 2 s d* p S /e 



dt 



12 



2tt 



exp 







6 



46 s 



2vr 

40 3 



exp 



6 



27(lnS) 2 _ 



27(lns) 



kT 
27rm! 



(S - l)A r 



27(ln S)' 



+ 2n 2 s v\ exp(0) 



kT 



{3-1) 



(52) 

(53) 
(54) 
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dt 27^(1115)2 



exp 







40 3 



27(lnS) 2 



kT 
27rmi 



Mi(5-1). (55) 



By eliminating dS/dt from Eq. (52) with the help of Eq. (53) and by applying a 

40 3 



finite difference formula in time, we get 



H t%— 1 



exp ^0 



12 V 2tt 



4e 3 



exp 



2/3 n 7 / 2 < 



27(lnS' i _ 1 )V 81 v^F5i 



27(ln5 s _ 1 )2_2 - jkT 



(ln^-i) 3 



5 



27rmi 



(56) 



Similarly, Eqs. (53)-(55) are discretized as, respectively, 



2/5 n 7 / 2 exp (e- 27(1 J i)2 ) n^Si-l)^ Hkf 



n\ 

51 V27T ^lllOi^i 
,2 



(In St. ^ 3 



Mi — Mi _./9u< /0 

-a \/ — exp 



H— 1 



+ 2v 



12Si p V 2vr 



40 3 



27rmi 



(57) 



27(lnSi_i) 2 



S? 



exp(0) 



£j «i— 1 



27V2t7 



(ln-Sj-i 



kT 
27rmi 



+ 4ttvi 



(5? 



2-K1TL 



(S t -1)M U . (59) 



By simplifying Eq. (56) with the help of Eqs. (57)-(59), we get a quadratic polyno 
mial of the form 

(bi -b 2 - h)n u 2 - n h + ni i _ 1 = , 



where 



h A A e 



b 2 = At 



40 3 
27(ln^_i) 2 

2(3 U 7 / 2 eX P ( Q ~ 27(lng.- 1 ) 2 



. Si - 1 kT 

b 3 = At * A n .\ 

d Sf ni V27rmi 



Solution of Eq. (60) yields two roots (denoted by superscripts (1,2)): 

1 



n 



(1,2) 

U 



2(6i - fe 2 - 6 3 



± 



1 - 4711^! (61 -b 2 -b 3 
2(6i - 62 - 63) 



(60) 
(61) 

(62) 
(63) 

(64) 
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For the first time step, the values of b\, 62 and 63 are obtained by applying the 
initial conditions: S(0) = So, n\ = n , and A no = A nQ . Since the n\ i must be real 
and finite, the following two conditions are imposed: 1 — 4ni < _ 1 (61 — 62 — 63) > and 
(61 — b 2 — 63) 7^ 0. Also, it has been assumed that the degradation of CNTs is an 
irreversible process, that is, the reformation of CNTs from the carbon cluster does 
not take place. Therefore, an additional condition of positivity, that is, > nx i _ 1 
is introduced while performing the time stepping. Along with the above constraints, 
the Hi history in a cell is calculated as follows: 

• If > rii t l and < n{ , then rii. = n^; 

• Else if n! 2 ' > ni. then ni. = n\ , 

• Otherwise the value of n\ remains the same as in the previous time step, that 
is, n u = n Xi _ x . 



Simplification of Eq. (57) results in the following equation: 
Si 2 + (ci + c 2 - Si-i)Si - ci = , 

where 



ci = Atn u A ni 



kT 
2irmi 



c 2 = At 



2(3 U 7 / 2 eX P ( - 27(lnf-i) 2 ) 
7^ n U 7T- o A3 



^ 1 (lnSV-x) 
Solution of Eq. (65) yields the following two roots: 

1 



Si 



? v ci + c 2 - ± ^\/ci + c 2 - Sf^ + 4ci 



(65) 

(66) 
(67) 

(68) 



For the first time step, c\ and C2 are calculated with the following conditions: 
from the above calculation, S(0) = S , and A no = A n0 . Realistically, the saturation 
ratio S cannot be negative or equal to one. Therefore, Si > yields C\ > 0. 



While solving for A n , the Eq. (59) is solved with the values of ri\ and S from the 

v4„o, Mi = M . The value of 



above calculations and the initial conditions A 



Mi was calculated by assuming n(d p , t) as a standard normal distribution function. 
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3.2 Incremental update of the CNT geometry 



At each time time step t = U, once the n\. is solved, we are in a position to compute 



the net electromechanical force (see Sec. 2.3) as 



f i = f i (E ,n li _ 1 ,e(ti-i)) 



(69) 



Subsequently, the orientation angle for each segment of each CNT is then obtained 



as (see Sec. 2.4) 



(70) 



and it is stored for future calculations. A critical angle, (9 C ), is generally employed 
with 9 C « 7i /A to 7r/2.5 for the present numerical simulations. For 9 < 9 C , the 
meaning of f z is the "longitudinal force" and the meaning of f x is the "lateral force" 



in the context of Eqs. (48) and (49). When 9 > 9 C , the meanings of f z and f x are 
interchanged. 



3.3 Computation of field emission current 

Once the updated tip angles and the electric field at the tip are obtained at a 
particular time step, we employ Eq. ([I]) to compute the current density contribution 
from each CNT tip, which can be rewritten as 

j, = —^\-—y (7i) 

with B = (1.4 x 1(T 6 ) x exp(9.8929 x $-V2) an d C = 6.5 x 10 7 taken from ref. 41 
. The device current (ij) from each computational cell with surface area AA at the 
anode at the present time step tj is obtained by summing up the current density 
over the number of CNTs in the cell, that is, 

I i = AA^2j i . (72) 
Fig. |5] shows the flow chart of the computational scheme discussed above. 
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At t — 0, in our model, the CNTs can be randomly oriented. This random 
distribution is parameterized in terms of the upper bound of the CNT tip deflection, 
which is given by Ax max = h/q, where h is the CNT length and q is a real number. 
In the numerical simulations which will be discussed next, the initial tip deflections 
can vary widely. The following values of the upper bound of the tip deflection have 
been considered: Aa; max = ho/ (5 + lOp), (p = 0, 1, 2, 9). The tip deflection Ax is 
randomized between zero and these upper bounds. Simulation for each initial input 
with a randomized distribution of tip deflections was run for a number of times and 
the maximum, minimum, and average values of the output current were obtained. 
In the first set, the simulations were run for a uniform height, radius and spacing 
of CNTs in the film. Subsequently, the height, the radius and the spacing were 
varied randomly within certain bounds, and their effects on the output current were 
analyzed. 

4 Results and discussions 

The CNT film under study in this work consists of randomly oriented multi-walled 
nanotubes (MWNTs). The film samples were grown on a stainless steel substrate. 
The film has a surface area of 1cm 2 and thickness of 10 — 14/xm. The anode consists 
of a 1.59mm thick copper plate with an area of 49.93mm 2 . The current- volt age 
history is measured over a range of DC bias voltages for a controlled gap between 
the cathode and the anode. In the experimental set-up, the device is placed within a 
vacuum chamber of a multi-stage pump. The gap (d) between the cathode substrate 
and the anode is controlled from outside by a micrometer. 

4.1 Degradation of the CNT thin films 

We assume that at t — 0, the film contains negligible amount of carbon cluster. 
To understand the phenomena of degradation and fragmentation of the CNTs, fol- 
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lowing three sets of input are considered: ni(0) = 100, 150, 500. The other initial 
conditions are set as 5(0) = 100, Afi(0) = 2.12 x 10~ 16 , A n (0) = 0, and T = 303K. 
Fig. [6] shows the three histories over a small time duration (160s) for the three 
cases of ni(0), respectively. For ui(0) = 100 and 150, the time histories indicate that 
the rate of decay is very slow, which in turn implies longer lifetime of the device. 
For ni(0) = 500, the time history indicates that the CNTs decay comparatively 
faster, but still insignificant for the first 34s, and then the cluster concentration 
becomes constant. It can be concluded from the above three cases that the rate of 
decay of CNTs is generally slow under operating conditions, which implies stable 
performance and longer lifetime of the device if this aspect is considered alone. 

Next, the effect of variation in the initial saturation ratio 5(0) on n\(t) history 
is studied. The value of ni(0) is set as 100, while other parameters are assumed to 
have identical value as considered previously. The following three initial conditions 
in 5(0) are considered: 5(0) = 50, 100, 150. Fig. [7] shows the n\(t) histories. It 
can be seen in this figure that for 5(0) = 100 (moderate value), the carbon cluster 
concentration first increases and then tends to a steady state. This was also observed 
in Fig. ^6J). For higher values of 5(0), ri\ increases exponentially over time. For 
5(0) = 50, a smaller value, the decay is not observed at all. This implies that a 
small value of 5(0) is favorable for longer lifetime of the cathode. However, a more 
detailed investigation on the physical mechanism of cluster formation and CNT 
fragmentation may be necessary, which is an open area of research. 

At t — 0, we assign random orientation angles (#(0) J ) to the CNT segments. 
For a cell containing 100 CNTs, Fig. [8] shows the terminal distribution of the CNT 
tip angles (at t = 160s corresponding to the ni(0) = 100 case discussed previously) 
compared to the initial distribution (at t = 0). The large fluctuations in the tip 
angles for many of the CNTs can be attributed to the significant electromechanical 
interactions. 
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4.2 Current-voltage characteristics 



In the present study, the quantum-mechanical treatment has not been explicitly 
carried out, and instead, the Fowler-Nordheim equation has been used to calculate 
the current density. In such a semi-empirical calculation, the work function <3> 42 
for the CNTs must be known accurately under a range of conditions for which the 
device-level simulations are being carried out. For CNTs, the field emission electrons 
originate from several excited energy states (non metallic electronic states) 43-44 . 
Therefore, the the work function for CNTs is usually not well identified and is more 
complicated to compute than for metals. Several methodologies for calculating the 
work function for CNTs have been proposed in literature. On the experimental 
side, Ultraviolet Photoelectron Spectroscopy (UPS) was used by Suzuki et al. 45 
to calculate the work function for SWNTs. They reported a work function value 
of 4.8 eV for SWNTs. By using UPS, Ago et al. 46 measured the work function 
for MWNTs as 4.3 eV. Fransen et al. 47 used the field emission electronic energy 
distribution (FEED) to investigate the work function for an individual MWNT that 
was mounted on a tungsten tip. Form their experiments, the work function was 
found to be 7.3±0.5 eV. Photoelectron emission (PEE) was used by Shiraishi et al. 48 
to measure the work function for SWNTs and MWNTs. They measured the work 
function for SWNTs to be 5.05 eV and for MWNTs to be 4.95 eV. Experimental 
estimates of work function for CNTs were carried out also by Sinitsyn et al. 49 . 
Two types were investigated by them: (i) 0.8-1.1 nm diameter SWNTs twisted into 
ropes of 10 nm diameter, and (ii) 10 nm diameter MWNTs twisted into 30-100 nm 
diameter ropes. The work functions for SWNTs and MWNTs were estimated to 
be 1.1 eV and 1.4 eV, respectively. Obraztsov et al. 50 reported the work function 
for MWNTs grown by CVD to be in the range 0.2-1.0 eV. These work function 
values are much smaller than the work function values of metals 3.6 — 5.4eV), 
silicon(^ 3.30 — 4.30eU), and graphite(^ 4.6 — 5.4eU). The calculated values of 
work function of CNTs by different techniques is summarized in Table I. The 
wide range of work functions in different studies indicates that there are possibly 
other important effects (such as electromechanical interactions and strain) which also 
depend on the method of sample preparation and different experimental techniques 
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used in those studies. In the present study, we have chosen $ = 2.2eV. 

The simulated current-voltage (I-V) characteristics of a film sample for a gap 
d = 34.7/im is compared with the experimental measurement in Fig. [9j The average 
height, the average radius and the average spacing between neighboring CNTs in the 
film sample are taken as h = 12/im, r = 2.75nm, and d\ = 1\im. The simulated I-V 
curve in Fig. [9] corresponds to the average of the computed current for the ten runs. 
This is the first and preliminary simulation of its kind based on a multiphysics 
based modeling approach and the present model predicts the I-V characteristics 
which is in close agreement with the experimental measurement. However, the above 
comparison indicates that there are some deviations near the threshold voltage of 
~ 500 — 600V, which needs to be looked at by improving the model as well as 
experimental materials and method. 



4.3 Field emission current history 

Next, we simulate the field emission current histories for the similar sample con- 
figuration as used previously, but for three different parametric variations: height, 
radius, and spacing. Current histories are shown for constant bias voltages of 440V, 
550V and 660V. 



4.3.1 Effects of uniform height, uniform radius and uniform spacing 

In this case, the values of height, radius, and the spacing between the neighboring 



CNTs are kept identical to the previous current-voltage calculation in Sec. 4.2 



Fig. 10 ^a), (b) and (c) show the current histories for three different bias voltages 
of 440V, 550V and 660V. In the subfigures, we plot the minimum, the maximum 
and the average currents over time as post-processed from a number of runs with 
randomized input distributions. At a bias voltage of 440V, the average current 
decreases from 1.36 x 10 _8 A to 1.25 x 10~ 8 A in steps. The maximum current varies 
between 1.86 x 10~ 8 A to 1.68 x 10~ 8 A, whereas the minimum current varies between 
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2.78 x 10- 9 A to 2.52 x 1(T 9 A Comparisons among the scales in the sub-figures 
indicate that there is an increase in the order of magnitude of current when the bias 
voltage is increased. The average current decreases from 1.25 x 10~ 5 A to 1.06 x 10~ 5 A 
in steps when the bias voltage is increased from 4401/ to 550V". At the bias voltage of 
660V, the average value of the current decreases from 1.26 x 10 _3 t4 to 1.02 x 10~ 3 A 
The increase in the order of magnitude in the current at higher bias voltage is due 
to the fact that the electrons are extracted with a larger force. However, at a higher 
bias voltage, the current is found to decay faster (see Fig. [T0[c)). 



4.3.2 Effects of non-uniform radius 

In this case, the uniform height and the uniform spacing between the neighboring 
CNTs are taken as ho = 12 nm and d\ = 2/im, respectively. Random distribution of 



radius is given with bounds 1.5— 4nm. The simulated results are shown in Fig. 11 At 
the bias voltage of 440V, the average current decreases from 1.37 x 10~ 8 v4 at t — Is 
to 1.23 x 10~ 8 A at t = 138s in steps and then the current stabilizes. The maximum 
current varies between 1.87 x 10~ S A to 1.72 x 10~ 8 A, whereas the minimum current 
varies between 2.53 x 10 _9 A to 2.52 x 10 _9 A The average current decreases from 
1.26 x 10~ 5 A to 1.08 x 10~ 5 A in steps when the bias voltage is increased from 440V 
to 550V. At a bias voltage of 660V, the average current decreases from 1.26 x 10~ 3 v4 
to 1.02 x 10~ 3 A As expected, a more fluctuation between the maximum and the 
minimum current have been observed here when compared to the case of uniform 
radius. 



4.3.3 Effects of non-uniform height 

In this case, the uniform radius and the uniform spacing between neighboring CNTs 
are taken as r = 2.75nm and d\ = 2fim, respectively. Random initial distribution 
of the height is given with bounds 10 — 14//m. The simulated results are shown in 



Fig. 12 At the bias voltage of 440V, the average current decreases from 1.79 x 10 6 A 
to 1.53 x 10 _6 A The maximum current varies between 6.33 x 10~ 6 v4 to 5.89 x 10~ 6 A 
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whereas the minimum current varies between 2.69 x 10~ 1C A to 4.18 x 10 _10 A The 
average current decreases from 0.495 x 10~ 3 v4 to 0.415 x 10~ 3 A in steps when the 
bias voltage is increased from 440V to 550V. At the bias voltage of 660V, the 
average current decreases from 0.0231A to 0.0178A The device response is found 
to be highly sensitive to the height distribution. 



4.3.4 Effects of non-uniform spacing between neighboring CNTs 

In this case, the uniform height and the uniform radius of the CNTs are taken as h = 
12 fim and r = 2.75nm, respectively. Random distribution of spacing d\ between 
the neighboring CNTs is given with bounds 1.5 — 2.5/zm. The simulated results are 



shown in Fig. 13 At the bias voltage of 440V, the average current decreases from 
1.37 x 10 _8 A to 1.26 x 10 _8 A The maximum current varies between 1.89 x 10~ 8 A 
to 1.76 x 10~ 8 A whereas the minimum current varies between 2.86 x 10~ 9 A to 
2.61 x 10 _9 A The average current decreases from 1.24 x 10~ 5 v4 to 1.08 x 10~ 5 t4 in 
steps when the bias voltage is increased from 440V to 550V. At the bias voltage of 
660V, the average current decreases from 1.266 x 10~ 3 v4 to 1.040 x 10~ 3 A There 
is a slight increase in the order of magnitude of current for non-uniform spacing. It 
can attributed to the reduction in screening effect at some emitting sites in the film 
where the spacing is large. 



5 Conclusions 



In this paper, we have developed a multiphysics based modelling approach to analyze 
the evolution of the CNT thin film. The developed approach has been applied to 
the simulation of the current- volt age characteristics at the device scale. First, a 
phenomenological model of degradation and fragmentation of the CNTs has been 
derived. From this model we obtain degraded state of CNTs in the film. This 
information, along with electromechanical force, is then employed to update the 
initially prescribed distribution of CNT geometries in a time incremental manner. 
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Finally, the device current is computed at each time step by using the semi-empirical 
Fowler-Nordheim equation and integration over the computational cell surfaces on 
the anode side. The model thus handles several important effects at the device 
scale, such as fragmentation of the CNTs, formation of the carbon clusters, and self- 
assembly of the system of CNTs during field emission. The consequence of these 
effects on the I-V characteristics is found to be important as clearly seen from the 
simulated results which are in close agreement with experiments. Parametric studies 
reported in the concluding part of this paper indicate that the effects of the height 
of the CNTs and the spacing between the CNTs on the current history is significant 
at the fast time scale. 

There are several other physical factors, such as the thermoelectric heating, 
interaction between the cathode substrate and the CNTs, time-dependent electronic 
properties of the CNTs and the clusters, ballistic transport etc., which may be 
important to consider while improving upon the model developed in the present 
paper. Effects of some of these factors have been discussed in the literature before 
in the context of isolated CNTs, but little is known at the system level. We note also 
that in the present model, the evolution mechanism is not fully coupled with the 
electromechanical forcing mechanism. The incorporation of the above factors and 
the full systematic coupling into the modelling framework developed here presents 
an appealing scope for future work. 
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Table 1: Summary of work function values for CNTs. 



Type of CNT 


$ (eV) 


Method 


SWNT 


4.8 


Ultraviolet photoelectron spectroscopy 45 


MWNT 


4.3 


Ultraviolet photoelectron spectroscopy 46 


MWNT 


7.3±0.5 


Field emission electronic energy distribution 


SWNT 


5.05 


Photoelectron emission 


MWNT 


4.95 


Photoelectron emission 48 


SWNT 


1.1 


Experiments 49 


MWNT 


1.4 


Experiments 49 


MWNT 


0.2-1.0 


Numerical approximation 50 



Figure 1: Schematic drawing of the CNT thin film for model idealization. 
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(a) (b) 

Figure 2: Schematic drawing showing (a) hexagonal arrangement of carbon atoms 
in CNT and (b) idealization of CNT as a one-dimensional elastic member. 



Figure 3: CNT array configuration. 
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Figure 4: Schematic description of neighboring CNT pair interaction for calculation 
of electrostatic force. 



Figure 5: Computational flow chart for calculating the device current. 
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Figure 6: Variation of carbon cluster concentration over time. Initial condition: 
S(0) = 100, T = 303fT, Mi(0) = 2.12 x lO" 16 , A n (0) = 0. 



Figure 7: Variation of carbon cluster concentration over time. Initial condition: 
ni(0) = 100m- 3 , T = 303K, Mi(0) = 2.12 x 10" 16 , A n (0) = 0. 
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Figure 8: Distribution of tip angles over the number of CNTs. 



Figure 9: Comparison of simulated current-voltage characteristics with experiments. 



(a) (b) (c) 

Figure 10: Simulated current histories for uniform radius, uniform height and uni- 
form spacing of CNTs at a bias voltage of (a) 440 V, (b) 550 V, and (c) 660 V. 
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(a) (b) (c) 

Figure 11: Simulated current histories for non-uniform radius of CNTs at a bias 
voltage of (a) 440 V, (b) 550 V, and (c) 660 V. 



(a) (b) (c) 

Figure 12: Simulated current histories for non-uniform height of CNTs at a bias 
voltage of (a) 440 V, (b) 550 V, and (c) 660 V. 



(a) (b) (c) 

Figure 13: Simulated current histories for non-uniform spacing between neighboring 
CNTs at a bias voltage of (a) 440 V, (b) 550 V, and (c) 660 V. 
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